MELLIN TRANSFORM AND IMAGE CHARGE METHOD FOR 
DIELECTRIC SPHERE IN AN ELECTROLYTE 
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Abstract. We revisit the image charge method for the Green's function problem of the Poisson- 
Boltzmann equation for a dielectric sphere immersed in ionic solutions. Using finite Mellin trans- 
formation, we represent the reaction potential due to a source charge inside the sphere in terms of 
one dimensional distribution of image charges. The image charges are generically composed of a 
point image at the Kelvin point and a line image extending from the Kelvin point to infinity with an 
oscillatory line charge strength. We further develop an efficient and accurate algorithm for discretiza- 
tion of the line image using Pade approximation and finite fraction expansion. Finally we illustrate 
the power of our method by applying it in a multiscale reaction-field Monte Carlo simulation of 
monovalent electrolytes. 
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1. Introduction. The method of image charges is a classical technique [38l [22] 
for electrostatic problems. Its most elementary application is the problem of a point 
charge in a spherical cavity inside a conducting media (or the reciprocal problem of a 
point charge outside a conducting sphere). In 1845, William Thomson (Lord Kelvin) 
[3"9"] noticed that the vanishing-potential boundary condition of conductors can be 
automatically satisfied on the sphere by putting an image point charge at the Kelvin 
point. A natural extension is the problem of a dielectric sphere inside a different 
I dielectric background, where a single point image charge no longer works. In 1883 

Carl Neumann [32] discovered that a point image at the Kelvin point together with 
a line image starting from the Kelvin point to infinity solves the boundary condition. 
This result has been independently re-derived by several authors, up to 1990s; see 
the reviews [321 HI] for more details. More recently, the problem of a spherical cavity 
inside an ionic solution has been studied 8 , 43 . It is again found that the image charge 
distribution consists of a point image and a line image. Computational method for 
the line charge density has been developed which works well in the asymptotic limits 
[8] [43] that kR is either large or small, where k is the inverse Debye length, while R 
is the radius of spherical cavity. The current work addresses the general case where 
kR is neither large nor small. 

The main advantage of image charge methods is to represent the effects of polar- 
ization charges in terms of point charges or line of charges, and therefore avoiding the 
task of numerically solving boundary value problems. This can substantially reduce 
the computational cost in Monte Carlo or molecular dynamics simulation of charged 
systems. Further reduction of computational cost can be achieved by discretizing 
line images using Gaussian quadratures, which effectively approximates a line charge 
by a few point charges. The total electrostatic energy of the system can then be 
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represented as Coulomb interaction between many point charges. 
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Fig. 1.1. Schematic illustrations of two multi-scale models of solute- solvent systems. The plus 
and minus signs represent the charges in the solute, the circles represent mobile ions in the solvent, 
and the angles represent the explicit water molecules, (a) The hybrid explicit/implicit solvent model. 
Inside the cavity, both water molecules and ions are treated explicitly. Outside the cavity, both 
bulk solvent and the ions are treated using the continuum theory; (b) The hybrid primitive/implicit 
solvent model. Inside the cavity, the water molecules are treated implicitly, while the ions are treated 
explicitly, i.e. the so-called primitive model. Outside the cavity, everything is treated implicitly. 

Most interfaces appearing in nature have irregular shapes, and the corresponding 
Green's functions can not be represented in terms of simple distributions of image 
charges. Hence one may rightfully argue that image charge methods have rather lim- 
ited applications. However, the spherical cavity does not have to be any dielectrics 
that is physically different from the exterior. Indeed, we believe that the most im- 
portant usage of image charge methods is in the multiscale reaction field modeling of 
Coulomb many body systems which arc of current interest in both simulations and 
continuum modelling [TTJ EH [H |I3 M ■ 

Due to the long range nature of Coulomb interaction, simulation of charged sys- 
tems is highly nontrivial. A proper treatment of boundary conditions is vital in order 
to obtain physically meaningful results. Periodic boundary conditions can remove 
artificial boundary effects in a self-consistent fashion and restore the translation sym- 
metry. Combined with Ewald summation method |12j . the cost of computing the 
total energy of the system with N particles scale as TV 3 / 2 . This can be further re- 
duced to the order of TV log N using a mesh-based algorithm such as the particle mesh 
Ewald or particle-particle particle-mesh Ewald lattice summation techniques. The 
periodic images are however unphysical and may produce artifacts that obscure the 
real physics. Besides, computational cost of Ewald-type summation method is still 
rather prohibitive for large systems in Monte Carlo simulations, which limits simula- 
tion of charged systems to rather small size. The development of non-Eward methods 
remains important topics; See Fukuda and Nakamura |17j for a recent review. 

An attractive alternative is to use the reaction field model, which is essentially a 
multi-scale strategy, schematically illustrated in Figurc fTTTl In this model, an artificial 
(spherical) cavity is introduced. All ions, together with possible mesoscopic objects 
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such as charged proteins and colloids, inside the cavity are treated explicitly, while 
the region outside the cavity is treated implicitly using appropriate continuum theory, 
such as linearized Poisson-Boltzmann theory. For any charge inside the cavity, then, 
we must solve the electrostatic boundary value problem, where the potential satisfies 
Poisson equation inside the cavity, and Poisson-Boltzmann equation outside. This 
problem can be efficiently solved using the image charge method developed in this 
work. We also note that there can be different levels of modeling inside the cavity. 
In the so-called hybrid implicit/explicit model [3S], both solvent molecules and ions 
are treated explicitly inside the cavity. By contrast, in the so-called hybrid primi- 
tive/implicit model, the solvent inside the cavity is modeled implicitly as a dielectric 
media, while the ions are treated at the level of primitive model. The primitive model 
is widely adopted in simulating soft matter systems |29j . 

The reaction field model augmented by our image charge methods can therefore 
provide a proper treatment of the electrostatic boundary conditions, but the computa- 
tion cost of the total energy scales as N 2 , where N here is the total number of particles 
including all ions as well as their images. Fortunately, this cost can be dramatically 
reduced using fast multipole based methods [201 111 HI HI] ■ The computational cost for 
the combined method generally scales as N. 

In this paper, we mainly focus on the image charge method for the Green's func- 
tion problem of these multi-scale reaction field models. The statistical mechanical 
foundation of these reaction field models will be discussed in a separate publication. 
In the remaining of this work, we shall first define the Kirkwood series for the Green's 
function (Sec. II) and then use the inverse Mellin transform to find the image charge 
representation (Sec. III). We further discretize the line image using method of Gauss 
quadrature and construct an efficient numerical scheme for the computation of the 
reaction potential. In Sec. IV, we compute the reaction potential using our method 
and quantify the errors. We also demonstrate the power of multi-scale modeling by 
Monte Carlo simulating a dilute symmetric electrolyte. 

2. Green's function of the Poisson-Boltzmann equation. As illustrated 
in Figure 12. 1[ let 17 £ R 3 be a spherical cavity with radius R, centered at the origin, 
with dielectric constant The volume outside the cavity is filled with electrolyte, 
characterized by a Debye length k~ 1 and dielectric constant e . For a source point 
r s inside the cavity, the relevant electrostatic Green's function satisfies the following 
equations: 

-£iV 2 $(r,r s ) =5(r-r s ), red, 

(2.1) 

— V 2 $(r, r s ) + K 2 $(r, r s ) = 0, otherwise, 

Throughout the paper we use light italic letter p to represent the magnitudes of a 
vector p. On the cavity boundary, the Green's function satisfies the following standard 
electrostatic interface conditions: 

$i = $ , £i3 r $i = £ 3 r $ , (2.2) 



1 Strictly speaking, grand canonical ensemble must be used in order to treat charge fluctuations 
properly. 

2 The parameters of the Poisson-Boltzmann theory however need to be determined self- 
consistently. Furthermore, extra short range interaction with the cavity boundary has also to be 
introduced to compensate the (unphysical) effects of the wall. Finally, to achieve a high efficiency, 
the radius of cavity should be chosen to be couple of the Debye length. 
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Poisson-Boltzmann equation 
Parameters: K, £ 
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Fig. 2.1. The geometry of the Green's function problem for a point charge at the source point 
r s ■ The reaction potential at a field point r is generated by a point image at the Kelvin point tk 
and a line image extends from to infinity. The line charge density is given by Theorem \3.1\ 



where $i and $ are the Green's function inside and outside the cavity respectively. 

The Green's function has an azimuthal symmetry, hence depends only on r, 6, 
where 9 is the angle between the source point r s and the field point r, see Fig. 12.11 
The potential inside the cavity, <£>i, is the superposition of the direct Coulomb potential 
'I'coui = l/(47T£i|r — r s |) and the reaction potential $ r f which is a harmonic function. 
Both potentials can be expanded in terms of spherical harmonics: 

li37T = E^r p "( C0S ^ ( 2 - 3a ) 

1 s 1 n— > 

oo 

$ rf (r,r s ) = Y^A n r n P n {cf»e), (2.3b) 

ra=0 

where r < (?* > ) is the smaller (larger) one between r s and r, while A n are constants to 
be determined, and P n ( - ) is Legendrc polynomial of order n. The potential outside 
the cavity can also be expanded in terms of Legendre polynomials: 

oo 

$ o = J2 B nkn(Kr)P n (cOs0), (2.3c) 
n=0 

where k n {-) is the modified spherical Hankel function [2] (also called the modified 
spherical Bcssel function of the second kind), defined by the following series, 



The coefficients A n ,B n in Eqs. ([2.3)1 can be found by solving the standard elec- 
trostatic boundary conditions Eq. (|2.2I) . We shall define the Kelvin point yk via 



v K = (R/r s ) 2 r s . (2.5) 
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For source charge inside the cavity, we always have 

r s < R< R 2 /r s =r K . (2.6) 

Using the orthogonality of the spherical harmonics, we obtain the following expansion 
of the reaction potential, widely known as Kirkwood series expansion |24j : 

$rf^E -s+rMn^Coosfl), (2.7) 

n— 

with the harmonic coefhcients given by 

MJu) = £ (" + 1 ) fc "H + 0) ^ (2 . 8 ) 
n enk n (u) — uk' n (u) 

where e = £\/s and u = kR. The Kirkwood series expansion has been used in the 
calculation of the reaction potential inside a spherical cavity in simulations |23| [6] . 
It converges slowly when a source charge approaches to the cavity surface, which 
prevents its wide application in dynamical and Monte Carlo simulations. 

3. Image charge representation and algorithm. In this section, we develop 
an image charge representation for the reaction potential using finite Mellin trans- 
formation. The Mellin transform method was used by Lindell and his collaborators 
28 , 33] for finding image charges of the Poisson equation. 



3.1. Finite Mellin transform. Let f(t) be a function defined in the interval 
[0, 1]. Its finite Mellin transform, F(s) is defined as [T5] 



l 

8-1. 



F(s)=M[f(t);s]= f(t)t s - L dt, (3.1) 
Jo 

where s is generically a complex variable. The original function f(t) can be expressed 
in terms of F(s) by the inverse Mellin transform 

f(t) = M- 1 [F(s);t}= F(s)r s ds. (3.2) 

J a — too 

Here the integration is carried out over a vertical line, Re(s) = a, in the complex 
plane. The dual functions {f(t),F(s)} form a finite Mellin transform pair. 

The finite Mellin transform is closely related to the one-side Laplace transform. 
Let t = e~ x , then A4[f(t); s] = £[f(e~ x ); s]. Likewise, the inverse transform can also 
be expressed in terms of the inverse Laplace transform. Two finite Mellin transform 
pairs shall be useful in our discussion below, 

{6(t-l), 1}, {t c , (s + c)- 1 }, (3.3a) 

where 6 is the Dirac delta function. The Mellin transformation of the second pair is 
defined for Re(s + c) > 0. 

We shall analytically continuate the integer variable n in Eq. (|2.8[) into the com- 
plex plane and let f(t) be its inverse finite Mellin transformation. Hence we have 

f(t) = M- 1 {M n {u);t), (3.4a) 
M n (u) = M(f; n)= [ f{t) t^dt. (3.4b) 
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Substituting Eq. (|3.4b[) into the Kirkwood series Eq. (|2.7[) . and using rpcr s = R 2 , the 
reaction potential can be re-expressed as 

$ ri = -r^r- [ dt r -^$- Y ^— T P„(cos0). (3.5) 

AnReJo i 2 ^ (r K /t) n+1 " l ' 1 ' 

Let us further define a vector x = yk /t. As t decreases from 1 to 0, the vector x runs 
from the Kelvin point rx to infinity along the radial direction. We shall see that this 
is precisely the loci of the image charge line. Let x be the magnitude of vector x, we 
have x = rx/t. Since we are only interested in the field point inside the cavity, we 
have r < R < rx /t = x, hence we can sum the series in Eq. (|3.5jl using the expansion 
Eq. (|2.3ap and obtain 

tl , ' r^^itsM. (3 . 6 ) 

47T£i J rK |r - x| 

The last integral represents the potential generated by one dimensional distribution 
of image charges along the radial direction, which starts from the Kelvin point yk and 
extends to infinity. The linear charge density is p(x) = R f{rn /x). The geometry 
is illustrated in Fig. 12.11 Therefore we arrive at our main result in this work. 

Theorem 3.1. Let {f(t),M n (u)} be a finite Mellin transform pair with t,n the 
conjugate variables, with M n (u) given in Eq. (|2.8[) . The Green's function problem 
(|2.ip has the following image charge representation: 

$ rf (r,r s ) = -^- f°° dx/^- (3.7) 

47T £i J rK \Y - X| 

where the line charge density is p{x) = f(rx/x)/R and x = (x/r s )Y s . 

The image charge representation for the reaction potential Eq. (|3.7[) is equivalent 
to the Kirkwood series representation, Eq. (|2.7[) . It is advantageous because the line 
integral can be efficiently discretized by Gauss quadrature. A few Gauss points can 
provide approximations with accuracy as high as desired [7] . Physically, this amounts 
to approximating line image by a few point images. 

3.2. Point image and line image. The following limit of the Bessel function 
can be established (using Eq. (|3.1T[) ) 

lim = (3 .8) 

n->oo n K n (U) 

Combining with Eq. (|2 .8[) we obtain the large n limit of the coefficients M n (u): 

7 = lim M n {u) = ^= £ -^-^. (3.9) 

n->co £ + 1 £■ + £ 



The coefficients Eq. (|2.8[) can therefore be decomposed into two parts 

M n {u) = 1 + 5M n {u), (3.10) 

where SM n (u) vanishes as n goes to infinity. 

The inverse Mellin transform of a constant 7 is a delta function 7<5(i — 1), while 
the inverse Mellin transform of the function SM n (u) is generally a continuous function 
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in the interval [0, 1]. The linear charge density therefore can be decomposed into (with 
t = r K /x): 

p{x) = R^M-^Mniu^t) 

= R-^Sft- 1) + R- 1 M- 1 {8M n {u)- : t) 

= ll S {x - r K ) + R^M- 1 (SM n {u); ^) . (3.11) 
r s \ x / 

The first term corresponds to a point image at the Kelvin point, and the second term 
corresponds to a continuous line image extending from the Kelvin point to infinity. 

The simplest limit is e — >• 0, where the exterior of the cavity becomes a conductor. 
In this limit, we easily see from Eq. (|2.8|l that M n (u) — > — 1 = 7, and SM n (u) = 
M„ (u) —7 = 0. The image charge distribution reduces to a single point image at the 
Kelvin point, a well known result. 

When the dielectric constants inside and outside the cavity are the same, i.e. 
s = l, the point image vanishes, but the line image persists. This is precisely the case 
of multi-scale reaction field model for electrolytes. 

3.3. Neumann's result revisited. Let us review the simple case where there 
is no screening ions outside the cavity, i.e. u = kR = 0. In 1883, Neumann [3"2"] 
found an exact expression for the reaction field, in terms of a point image charge at 
the Kelvin point and a line image: 

a _ iR/rs , is [°°, (rWx) 1/(E+1) ,„ 1<rt 
* rf - 4n £i \r r K \ + (e + 1)R J ric 47r £i |r - x| ' (6AZ) 

where 7= (e — l)/(e + 1). This formula has been re-derived independently by various 
authors in different fields of applications [44] [13l [TUl EH [27l [34]; also see Lindell's 
review [36] for the summary of history. 

The Neumann's integral expression can be easily derived from the general result 
Eq. (|3.7I) . For u = 0, the harmonic coefficients Eq. (|2.8[) reduce to, 



M„(0)=7+ Z't+l, 
n + (e + 1) 1 

Its inverse Mellin transform can be exactly calculated using Eq. ([3.3 



f(t) = 1 6(t-l) + j^t 1 ^+ 1 l (3.13) 

Substituting this back into Eq. (|3 . T[) we find Neumann's result Eq. (|3.12l) . 

The application of the Neumann's result in molecular dynamics can be found in 
two recent papers [13 [55] . The reciprocal problem of a source charge placed outside 
of a sphere was also applied in Monte Carlo simulations of colloidal systems [T5] . 
Historically, what is widely used in computer simulations of biological systems is only 
single image charge approximations, see works by Friedman |16| and Abagyan and 
Totrov [T] . These methods are of the first or second order accuracy in the dielectric 
ratio e = £i/e - They fail to be accurate if the ratio e is not small. 

3.4. Limits of large cavity and small cavity. For u = kR ^ 0, the inverse 
Mellin transform of M s (u) can not be exactly calculated. One way to proceed is to 
use the following asymptotic approximation [43] . 



uk' n (u) 
k n (u) 



1 
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which gives the correct leading order asymptotics both for u — > and for u — > +00. 
Substituting it back into Eq. (|2.8[) leads to 

1 -2e(l + a)/(l+g) 
M " (U) " 7+ (l +g )n+(l + fi) ' (3 ' 14) 

where {t = u 2 / (l + u). The inverse Mellin transform of this can be easily found. This 
approximation works well for the case of small u (small cavity) and large u (large 
cavity) . 

3.5. Large n asymptotics of M n (u). Let us first look at the large n limit of the 
Bessel function k n (u). For sufficiently large n, the largest term in the sum Eq. (|2.4I) 
is given by I = n. We can therefore rewrite the summation as 

_ ffe -"(2n)! ^ n!(2n~p)! 
MU) " (2u)™+%! ^ (2n)!(n - p)\p\ [2u) ' ^ 15j 

where we have defined p = n — I. We can expand the function being summed into 
asymptotic series in terms of 1/n and extend the upper limit of summation from n to 
00. The resulting summation then can be calculated order by order in 1/n. This can 
be conveniently done using Wolfram Mathematica. For example, up to order of n~ 2 , 
we have 

2-"-6 u -«-ir(2n + 1) Su 2 ( u 4 -4m 2 



k n (u) = - - \ 32 + ^ + (n 3 ) (3.16) 

1 (n + 1) \ n / 

It then follows that 

, x (e-l)e e(u 2 e + u 2 + e-l) . „. . 

M " w = 7 ,12 - 2TTTV3 " + ° " r ' 3 ' 18 

n(e + \y n A [€ + 1) J 

Generically, therefore, the leading order term of M n {u) scales as 1/n. In the most 
interesting case e = 1, however, this term vanishes and we have 

u 2 

M n (u) = -^ + 0(n- 3 ). (3.19) 

Now consider the limit where the source charge approaches the cavity boundary, we 
have r s — > R, tk — > R- The reaction potential acting on the source charge then (see 
Eq. ([277]) ) becomes 

$ rf (R, R) -> V M n («) < 00. (3.20) 

Therefore the line image strength must vanish at the Kelvin point. 
3.6. The general case. 
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3.6.1. Pade approximations to harmonic coefficients. In the multi-scale 
reaction field model, we typically have e«1 and u = kR of order of unity. All methods 
discussed above fail in this case. To obtain accurate approximation, we approximate 
the harmonic coefficients M n {u) by a rational function of n (Pade approximation): 



M(P; n) = 7 + V / V , (3.21) 




with constants oij £ K and j3j > for all j. Clearly Eq. f|3 . 14[) is a special case 
of Eq. p.21j) with P = 1. The second term on the right hand side is the Pade 
approximation of the function SM n (u) of order (P — 1, P) , and the rational polynomial 
preserves the asymptotics, 

M n (u) =7 + 0(l/n), as n -> oo. (3.22) 

The coefficients {a^, /3j} are solved by the nonlinear least square method. For given u 
and dielectric ratio e, the constants in the expansion can be simply determined by a 
minimization of the total L2 error of the first N + 1 terms of the harmonic coefficients, 

N 



min V \M n (u) -M(P;n)] . (3.23) 



The Newton iteration scheme or other iteration algorithms can be applied to solve 
this nonlinear optimization problem. 

To demonstrate the quality of Pade approximation, we list in Table lBTTI thc relative 
errors for the cases of P = 2 and 3 with various parameters u and e: 

N 2 / N 

E 2 = Y, [M n (u) - M(P; n)] / £ M n {uf , 

n=0 I n=0 

where M(P; n) is the numerical solution of the nonlinear least square problem (|3.23l) . 
It should be pointed out that the minimization solution depends on the initials and 
may be not a global minimization. Nonetheless, the results in Table l3~T1 clearly show 
the remarkable precision of the Pade approximation. For example, the worst case for 
P = 3 has a relative error ~ 0.02%. In comparison, the asymptotic method Eq. (|3.14l) 
yields much larger error, and completely breaks down for e — > 1. 

3.6.2. Image expressions through inverse Mellin transforms. After the 
best- fitting coefficients {ctj,[3j} are determined, we can further re-express the Pade 
approximation Eq. (|3 . 2 1 1) using the partial fraction expansion: 

P c 

M(P;n)= 7 + ^— V- ( 3 - 24 ) 

Using Eq. (|3.3p , the inverse Mellin transform f(t) = M^ 1 [M(P; s);t], we easily find 
the image charge distribution: 

p(x) = -^S(x-rK) + X(r K /x), (3.25) 
r. 
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Table 3.1 

Relative errors of Fade approximations to the harmonic coefficients for different parameters u 
and e. In the minimization, N = 50 is taken. The asymptotic solution Eq. \3. 14]) ctlso given 

for comparison. 



£ = 0.02 £ = 0.1 £ = 0.5 £ = 1 

P = 2 

w = 0.5 4.10E-6 2.03E-5 1.10E-4 6.01E-4 
u = 2 3.14E-5 1.64E-4 1.15E-3 5.07E-3 
u = 10 6.08E-5 3.18E-4 2.08E-3 5.74E-3 

P=~3 

u = 0.5 3.56E-9 2.17E-8 4.25E-7 1.16E-5 
u = 2 9.55E-7 4.92E-6 3.31E-5 1.41E-4 
u — 10 2.22E-6 1.15E-5 7.37E-5 1.99E-4 

Asymptotics Eq. (|3.14|) 

u = 0.5 1.12E-3 5.84E-3 4.28E-2 0.42 

u = 2 1.67E-3 9.34E-3 9.05E-2 0.58 

u = 10 8.89E-3 4.66E-2 0.32 0.95 



where 

i p 

3=1 

with t — rji/x represents a line image density from the Kelvin image point to the 
infinity along the radial direction. Alternatively, the reaction potential is, 

<Mr,r s ) = ^ : + — r dx KrK,X } . (3.27) 

rU s > 4irRei\r-r K \ inei J rK |r - x| y J 

Figure [37X1 illustrates the profiles of line image strengths A(rx/x) for u = 2 and for 
£ varying from 0.02 to 1 calculated with P = 3. Note that the line image is always 
oscillatory, which implies that some of the parameters Zj are complex numbers. Note 
also that both amplitude and period increase with the dielectric ratio e. 

3.6.3. Discretization of the line image. Computation of the reaction poten- 
tial using the continuous linear charge density Eq. ()3.26j) is still expensive. Therefore 
we further discretize the line integral (|3.27[) using numerical quadrature. This amounts 
to approximating the linear image using multiple point image charges. An efficient 
discretization scheme uses fewer effective point images and hence significantly reduces 
the computational cost in simulations. 

We consider the case of P = 3. Extension of the algorithm to other values of 
P is straightforward. Since M(P; n) must be real and the line charge is oscillatory, 
the parameters {Zi, Z 2 , Z 3 } in Eq. Q3.24[) are generally composed of a positive real 
number and a complex conjugate pair. The same clearly also holds for the set of 
parameters {Ci, C2, C3}. Let Z\,C\ be real and 

Z 2 , 3 =p±iq, C 2 , 3 = a±ib. (3.28) 

The line image strength in Eq. f|3 . 26[) can be expressed as 

\{t)R= Ci t Zl +2t p [a cos(g log t) + b sin(g log t)} , (3.29) 
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Fig. 3.1. Profiles of line image strengths as a function of x/rx for parameters u = 2 and 
different dielectric ratios £ calculated with P = 3. Note that the linear charge density is always 
oscillatory in radius. 



where t = tk /x. Note that the second term is oscillatory, in agreement with Fig. 13.11 
Approximating an oscillating integral is tricky. We divide the integral in Eq. (|3.27l) 
into two parts: 



r-x 



\{r K /x) 



r-x 



dx 



\{r K /x) 



dx. 



r K T 



r-x 



(3.30) 



with T a positive number sufficiently larger than unity. The first integral in Eq. (|3 . 30[) 
can be transformed into an integral over the interval v 6 [— 1, 1] through the linear 
variable transformation: 

x(v) = r K (2 + e) T /(l + e-v) T 
with e = 2/ (T 1 / 1 " — 1) and r a positive constant. Further defining a function Q(v) via 

tx{v) 



Q(v) = \(r K /x(v)) 



1 + e-v' 



we can easily show that 



" T ^4 dx .r 9w dv . 

|r-x| J_ x |r-x(u)| 



(3.31) 



(3.32) 



The resulting integral over v can be then integrated using the classical M-point Gauss- 
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Legcndre quadrature, which discretizes the line image into several point image charges 

"""'^E^ (3-33) 



where 



r - x ' r - r ri 

1 1 m— 1 



— -^(^m)i Qm — Q{,^m) 



and {s m , u> m , m — 1, • ■ • , M} are the Gauss quadrature points and weights. 

For the second integral in the RHS of Eq. p.30p . the variable x is much larger 
than r, hence the integrand can be expanded in terms of the ratio r/x: 



r — x i — ' x 

1 1 n=0 



Therefore 



I 



X(r K /x) 



dx = -J2lir n P n {cos9) 7 (3.34) 



where 



r ^T r — x R 



Jr K T x 



The first few terms can be explicitly worked out using Eq. (|3.29[) : 



C\T 


-Zi 

h 


2 (a P 


Z 


l 




±\ 


\C{T~ 


-Zi-l 


r K 


{ Zi 


+ 1 




\G X T~ 


-Zi-2 


^\ 


[ Zi 


+ 2 



TP+![(p+l) 2 + q 2 

^ [o(p + 2) - bq\ c 



TP+ 2 [(p + 2) 2 + g 2 ] 

(3.35) 

For sufficiently large T, the zeroth order correction Iq already gives an approximation 
with high accuracy. 

In summary, the reaction potential Eq. (|3.27|) is approximated by M + 1 image 
point charges plus a few correction terms, 

^-^Xv^\ + ^R^ Iirn Pnicos 91 (3 - 36) 

m=0 1 1 1=0 

where the term of m = represents the Kelvin image charge, with q$ = yrx/R, *o = 
tk- As pointed out previously, the Kelvin image vanishes when e = 1. 

4. Numerical results. In this section, we compute the reaction potential using 
our method and quantify the errors. We also demonstrate the power of multi-scale 
modeling by Monte Carlo simulating a dilute symmetric electrolyte. 
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4.1. Accuracy of discretization scheme. We first test the accuracy of our 
discretization scheme. We focus on the most difficult case e = 1 where the line image 
is strongly oscillatory, and the point image vanishes. It is also the case appears in 
the multiscale modeling of electrolytes. We take u = 5 and R = 1, and calculate 
the relative error of the self energy of a unit charge, $ r f (r s , r s )/2, as a function of 
the source point radius r s , by comparing with 201-terms-truncation of the Kirkwood 
series. We set r = 5 and T — 4. We find that the accuracy of numerical quadrature 
only weakly depends on T and r. 

We use 4,5,6, and 15 Gauss quadrature points for the integration on the finite 
interval \tk , trT] and the two different corrections for the integration on interval 
[rj(T, oo): L = and 1. Note that L = means the second term of Eq. (|3.36j) is a 
constant correction. The results are shown in Figure 14.11 It is seen that the image 
charge approximation is very accurate even with only 4 image points, with the overall 
relative errors remaining less than 1%. The accuracy however does decrease when r s 
approaches to the boundary. 

To determine the asymptotics of image charge approximation for the self energy 
near the boundary, we also compare it with the direct truncation of the Kirkwood 
series, in the range r s /R G [0.9,0.99]. These results, shown in Figure |4~2| clearly 
demonstrate that image charge approximation converges much faster than the Kirk- 
wood series near the boundary. In particular, image charge approximation with 4 
image charges is uniformly better than Kirkwood series with 20 terms. For example, 
the error of Kirkwood series with 20 terms is larger than 1% for r s /R = 0.97, and 
increases to 4% for r s /R = 0.99, while the error of the image charge approximation 
remains less than 1%. 

4.2. Application in Monte Carlo simulations of electrolytes. To illustrate 
the utility of our image charge method, we apply it in a reaction field Monte Carlo 
simulation of electrolyte. We artificially introduce a spherical cavity with radius R 
in the electrolyte, and model all charges in the cavity using the primitive model 
and simulate them using Monte Carlo method. All ions outside the cavity, together 
with the solvent, are treated implicitly using linearized Poisson-Boltzmann theory, 
characterized by two parameters: the dielectric constant e Q and the inverse Debye 
length k. For any ion inside the cavity, the effects of the electrolyte outside the cavity 
is to introduce a reaction potential, which can be calculated using our image charge 
method. It is important to note that the inverse Debye length k characterizing the 
medium outside the cavity shall be self-consistently determined by the simulation of 
the ions inside the cavity. 

We run canonical Monte Carlo simulations of the primitive model using the stan- 
dard Metropolis criterion [3TJEI] for particle displacements. The ions are modeled by 
hard spheres with diameter D = 3. 75 A and with a point charge of valence Zi = ±1 
at its center. The ions are mobile in the solvent medium with dielectric permittivity 
Ej = 80. The effective Hamiltonian of the system is given by 



N 




(4.1) 



The self energy 



is given by 




n < R, 
n > R, 



(4.2) 
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Fig. 4.1. Relative errors of the self energy as a function of source position r s for different 
numbers of images (A, 5, 6, 15). (a) L = 0; (b) L = 1. 



where Ib — e 2 / (4irEoEikBT) is the Bjerrum length and /3 = l/(ksT) is the Boltzmann 
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Fig. 4.2. Relative errors of the self energy as a function of source position r s 6 [0.9,0.99] for 
different numbers of images (4 ~ 6) with constant correction L = and different numbers of poles 
with multipole expansion. 



factor. We take Ib = 7.14A for the water permittivity at room temperature. The 
infinite potential is due to the presence of the hard wall at r = R + D/2. The pairwise 
interaction energy [Zy is given by 



f3Uij = IsZiZj 



- — - +47re i $ rf (r i ,r i ) 



(4.3) 



As is well known, the reaction potential is symmetric under permutation of two vari- 
ables: $ rf (ri,rj) = $rf(r,-,ri). 

We first calculate the density distribution by taking a system with R = 100A and 
the salt concentration is po = 8mM which corresponds to a Debye length k" 1 = 34 A 
We calculate both the cases with and without the reaction field. In the former case, 
we use different numbers of image charges varying from 4 to 6, with the parameters 
for numerical quadratures the same as those in the upper panel of Figure 14.11 For 
each setting, we run 6 x 10 8 MC cycles for each particle to obtain samples for the 
statistics of particle number in each spherical shell with thickness 2A. 

In Figure I4TB1 we show the radial distribution function of anions, i.e., the nor- 
malized density, p/po- The distribution of cations is similar. Evidently, without 
accounting for the reaction field, the density is higher near the center of cavity and 
lower near the cavity boundary. The difference in the density is up to 3 percent. 
When the reaction field is taken into account, the particle density shows variation less 
than 0.5% in most regions inside the cavity. There is substantial deviation of density 
near the cavity boundary. This is due to the presence of the artificial hard wall on 
the boundary, sec Eq. (|4.2[) . This problem can be cured by adding to the self energy a 
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Fig. 4.3. Comparison of radial density ratios of anions in the spherical cavity for approximating 
the reaction field with different numbers of images, and without reaction field. Cavity radius is 
R = lOOA; the salt concentration is po = &mM , corresponding to a Debye length k~ 1 = 34/i. 



short ranged correction, or by introducing a buffer zone of a certain thichness (~ D). 
The latter approach has been discussed in literature [3J 2U HI]- We shall present 
statistical mechanical discussion of the former approach in a separate publication. 

Finally, in Figure 14.41 we plot the results for three different bulk ionic concen- 
trations, where p = 8, 16 and 24mM respectively. These results again show the 
necessity of treating the reaction potential, as well as the accuracy and efficiency of 
our image charge methods. It is also interesting to note that as the Debye length 
decreases the density becomes flatter, suggesting that the reaction field Monte Carlo 
models becomes more precise. For the case of salt concentration 24mM, for example, 
the Debye length is approximately half of the cavity radius, while the density variation 
becomes less than 0.1% (excluding the thin region affected by the hard wall artifacts). 

5. Conclusions. In summary, we have developed an image method for charges 
inside a spherical cavity that is immersed in an ionic solution. Our method is useful 
for multiscale reaction field models of electrolytes and other more complicated charged 
systems. We derive an analytic expression for the reaction potential in terms of one 
dimensional image charge distribution, and discuss a highly accurate and efficient 
algorithm for discretizing the image line charge. We also apply our method to a 
reaction field Monte Carlo simulations of 1 : 1 electrolytes, and demonstrate the 
accuracy and efficiency of the new algorithm. 

Simulation of charged systems is computationally expensive, therefore is always 
limited to small system size. In a physical system of such size, the total charge 
may fluctuate away from zero to a noticeable extent. These fluctuations can not 
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be taken into account in canonical ensemble simulations. In another word, grand 
canonical ensemble must be used to capture the charge fluctuations of small systems. 
Ewald-summation method, which is so far the most popular simulation methods for 
charged systems, are based on periodic boundary conditions, and are difficult to be 
incorporated with grand canonical ensemble. By strong contrast, reaction- field type 
of models, besides being more intuitive, can be easily adapted to a grand canonical 
Monte Carlo simulation. This shall be the topic of a separate publication. 
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